 #This file estimate EV
# Imports
from coopr.pyomo import *
from coopr.opt import SolverFactory
from ReferenceModel import model
import numpy

# Solve WS for given number of sample realizations
numSamples=20
numX = 5

optVal  = numpy.array ([0 for i in range(numSamples)])
optSoln = numpy.array([[0 for i in range(numSamples)] for j in range(numX)])

opt = SolverFactory('gurobi')

for i in range(numSamples):
    datafile = './scenariodata/Scenario' + str(i+1) + '.dat'

    instance = model.create(datafile)

    results = opt.solve(instance)
    print "Solve" + str(i) + "th instance"
    instance.load(results)
    optVal[i] = value(instance.TotalProfit)

    for j in range(numX):
        optSoln[j][i] = value(instance.X[j+1])

# Calculate point est / interval est of objective value, point estimate of solution
z_alpha = 1.96

WS     = optVal[:].mean()
WS_var = optVal[:].var() * numSamples/(numSamples-1)
WS_halfwidth = z_alpha * sqrt(WS_var/numSamples)

WS_CI_lo = WS - WS_halfwidth
WS_CI_hi = WS + WS_halfwidth

X_WS = [0 for j in range(numX)]

for j in range(numX):
    X_WS[j] = optSoln[j][:].mean()

print WS
print WS_CI_lo, WS_CI_hi
print X_WS


